######################################################################################################################################
####### The code below pertains to tables (1-2) and Figure (4), constituting the main analysis presented in the paper ################
######################################################################################################################################

# Load the required packages
rm(list = ls())
library(plyr)
library(dplyr)
library(ggplot2)
library(data.table)
library(caTools)
library(haven)
library(readr)
library(dummies)
library(lubridate)
library(fixest)
library(readstata13)
library(robustHD)
library(AER)
library(stdidx)
library(jtools)
library(MASS)
library(mosaic)
library(sensemakr)
library(stdidx)
library(coefplot)
library(collapse)

#Loading the two main datasets: 
# data_2rounds: Contains data on the 2 rounds of the election, excluding the pre-electoral period 
# data_3rounds: Contains data on the 2 rounds plus the pre-electoral period as a "third round" 

setwd("")
data <- read_dta("data_2rounds.dta")
data3 <- read_dta("data_3rounds.dta")

#Create dictionary to present the names of the variables in the tables
dict=c("mb_dummy"="MB Running", "pct_urban"="Urban (%)", "pct_emp"="Employment (%)", 
       "sdeduscore"="Education (sd)","eduscore2"="Education2 (sd)", "gov_employed"="Gov. Employment (%)", "logprotest"="Protest (log)", 
       "highprofile"="High Profile", "logregistered"="Registered (log)", "female_pct"="Female (%)", "round2"="Runoff", 
       "votebuyingd"="Vote-buying", "fmean.female_pct"="Female", "fmean.logregistered"="Registered (log)", 
       "fmean.pct_urban"="Urban (%)", "fmean.pct_emp"="Employment (%)", "fmean.sdeduscore"="Education (sd)", 
       "fmean.gov_employed"="Gov. Employment (%)", "newndp"="New NDP", "split_ndp"="NDP Dissidents", 
       "incumbent"="Incumbents", "nocand_competing"="Candidates No.", "death_ndpaffil"="NDP Death", "client_index"="Clientelism", 
       "watani_1984" = "NDP 1984 (%)", "wafd_1984" = "MB Coalition 1984 (%)", "fmean.incumbent"="Incumbents", "fmean.nocand_competing"="Candidates No.", 
       "fmean.newndp"="New NDP", "fmean.split_ndp"="NDP Dissidents", "fmean.logprotest"="Protest (log)")

####################### Main Models for the Paper #########################################

######## Table (1)
main1 <- fenegbin(allintimid ~ newndp +  split_ndp + mb_dummy + pct_emp + sdeduscore + pct_urban   + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing + round2 | governorate, data)

main2 <- fenegbin(allintimid ~ newndp +  split_ndp + mb_dummy +  pct_emp + sdeduscore +pct_urban  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==1))            

main3 <- fenegbin(allintimid ~ newndp +  split_ndp + mb_dummy + pct_emp + sdeduscore +pct_urban  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==2))     

main4 <- fenegbin(allintimid ~ newndp +  split_ndp + mb_dummy + pct_emp + sdeduscore +pct_urban   + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing+ round2+round1 | governorate, data3)            


etable(main1, main2, main3, main4,  se = c("cluster"), cluster=~governorate, digits = 3, tex=T, sdBelow = T, 
       signifCode = c("***"=0.001, "**"=0.01, "*"=0.05, "+"=0.10), dict=dict, interaction.combine = "*")



######### Figure (4)

# Panel (a): Violence against voters

#Estimation

mainv1 <- fenegbin(voter_total ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing + round2 | governorate, data)

mainv2 <- fenegbin(voter_total ~ newndp +  split_ndp +mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==1))            

mainv3 <- fenegbin(voter_total ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==2)) 

mainv4 <- fenegbin(voter_total ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing+ round2+round1 | governorate, data3)          

#Coefficient plot

fixest::coefplot(list(mainv1,mainv2, mainv3, mainv4), 
                 keep=c("MB Running", "New NDP", "NDP Dissidents", "Urban", "Employment", "Education", "Protest"), 
                 dict=c("reset", "mb_dummy"="MB Running", "pct_urban"="Urban (%)", "pct_emp"="Employment (%)", 
                        "sdeduscore"="Education (sd)","newndp"="New NDP", "split_ndp"="NDP Dissidents", "logprotest"="Protest"), 
                 lab.fit="multi", horiz=T, main="") 

legend("bottomright", col = 1:4, lwd = 1, pch=c(16,17,15,13),
       legend = c("2 Rounds", "First", "Second", "3 Rounds"), title = "Model")



# Panel (b): Violence by state actors

#Estimation

mains1 <- fenegbin(state_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing + round2 | governorate, data)

mains2 <- fenegbin(state_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==1))            

mains3 <- fenegbin(state_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==2))            


mains4 <- fenegbin(state_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing+ round2+round1 | governorate, data3)          


#Coefficient Plot

fixest::coefplot(list(mains1,mains2, mains3, mains4), 
                 keep=c("MB Running", "New NDP", "NDP Dissidents", "Urban", "Employment", "Education", "Protest"), 
                 dict=c("reset", "mb_dummy"="MB Running", "pct_urban"="Urban (%)", "pct_emp"="Employment (%)", 
                        "sdeduscore"="Education (sd)","newndp"="New NDP", "split_ndp"="NDP Dissidents", "logprotest"="Protest"), 
                 lab.fit="multi", horiz=T, main="") 

legend("bottomright", col = 1:4, lwd = 1, pch=c(16,17,15,13),
       legend = c("2 Rounds", "First", "Second", "3 Rounds"), title = "Model")



# Panel (c): Violence by non-state actors

#Estimation 

mainns1 <- fenegbin(nonstate_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing + round2 | governorate, data)

mainns2 <- fenegbin(nonstate_intimidation ~ newndp +  split_ndp + mb_dummy + pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==1))            

mainns3 <- fenegbin(nonstate_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing | governorate, subset(data, data$round==2))  

mainns4 <- fenegbin(nonstate_intimidation ~ newndp +  split_ndp + mb_dummy +  pct_urban + pct_emp + sdeduscore  + logprotest + logregistered+female_pct
                  + incumbent + nocand_competing+ round2+round1 | governorate, data3)          


#Coefficient plot 

fixest::coefplot(list(mainns1,mainns2, mainns3, mainns4), 
                 keep=c("MB Running", "New NDP", "NDP Dissidents", "Urban", "Employment", "Education", "Protest"), 
                 dict=c("reset", "mb_dummy"="MB Running", "pct_urban"="Urban (%)", "pct_emp"="Employment (%)", 
                        "sdeduscore"="Education (sd)","newndp"="New NDP", "split_ndp"="NDP Dissidents", "logprotest"="Protest"), 
                 lab.fit="multi", horiz=T, main="") 

legend("bottomright", col = 1:4, lwd = 1, pch=c(16,17,15,13),
       legend = c("2 Rounds", "First", "Second", "3 Rounds"), title = "Model")


############# Table (2)

#This analysis requires: 

# (1) Aggregating the main variables for the two rounds

data2 <- collap(data, allintimid+mb_dummy + pct_urban+pct_emp + sdeduscore  + logprotest +
                  logregistered+female_pct+round2+incumbent+nocand_competing+newndp+split_ndp ~ governorate+district_name, list(fsum, fmean))

# (2) Matching the violence data to clientelism data
clientdata <- read_dta("clientelismdata.dta")
data2 <- merge(data2, clientdata, by=c("district_name"))

# (3) Adjusting the definitions of two variables
data2$mb_dummy <- as.numeric(data2$fmean.mb_dummy>0)
data2$votebuyingd <- as.numeric(data2$votebuying>0)


#Finally, we estimate the data for Table (2)

m1 <- fenegbin(fsum.allintimid ~ votebuyingd+ fmean.female_pct + fmean.logregistered +fmean.incumbent+fmean.nocand_competing| governorate, 
               data2)

m2 <- fenegbin(fsum.allintimid ~ votebuyingd + mb_dummy+fmean.newndp+fmean.split_ndp+
                 fmean.logprotest+ fmean.female_pct + fmean.logregistered +fmean.incumbent+fmean.nocand_competing| governorate, 
               data2)

m3 <- fenegbin(fsum.allintimid ~ votebuyingd + mb_dummy+fmean.newndp+fmean.split_ndp++ fmean.pct_urban + fmean.pct_emp + fmean.sdeduscore+
                 fmean.logprotest+ fmean.female_pct + fmean.logregistered +fmean.incumbent+fmean.nocand_competing| governorate, 
               data2)

m4 <- fenegbin(fsum.allintimid ~ votebuyingd + mb_dummy+fmean.newndp+fmean.split_ndp++ fmean.pct_urban + fmean.pct_emp + fmean.sdeduscore+
                 fmean.logprotest+ fmean.female_pct + fmean.logregistered +fmean.incumbent+fmean.nocand_competing, 
               data2)

etable(m1, m2, m3, m4,  se = c("cluster"), cluster=~governorate,  digits = 3, tex=T, sdBelow = T, 
       signifCode = c("***"=0.001, "**"=0.01, "*"=0.05, "+"=0.10), interaction.combine = "*",  dict=dict)

